A fast algorithm for approximating the ground state energy on a quantum computer 



Anargyros Papageorgiou/'B Iasonas Petras, 1 '!!) Joseph F. Traub/'p and Chi Zhang 1 '!! 

1 Department of Computer Science, Columbia University, New York, NY 10027, USA. 

(Dated: October 12, 2010) 

Estimating the ground state energy of a multiparticle system with relative error e using deter- 
ministic classical algorithms has cost that grows exponentially with the number of particles. The 
problem depends on a number of state variables d that is proportional to the number of particles and 
suffers from the curse of dimensionality. Quantum computers can vanquish this curse. In particular, 
we study a ground state eigenvalue problem and exhibit a quantum algorithm that achieves relative 
error e using a number of qubits C'dloge -1 with total cost (number of queries plus other quantum 
operations) Cde~^ +S \ where 8 > is arbitrarily small and C and C' are independent of d and e. 
Thus, the number of qubits and the total cost are linear in the number of particles. 
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Introduction. A difficult and challenging problem in 
modern science is to accurately compute properties of 
physical and chemical systems. One of the difficulties 
in carrying out precise calculations arises from the com- 
putational demands the Schrodinger equation presents. 
The computational resources needed to obtain accurate 
solutions appear to be exponential in the size of the phys- 
ical system. As a result these problems are considered in- 
tractable on a classical computer for systems that are not 
trivial in size. For an overview of the numerical methods 
used for the solution of such problems see [l], and the 
references therein. 

So far there have been mixed results about the po- 
tential power of quantum computers relative to that of 
classical computers. For some problems, such as factor- 
ing large numbers, quantum computers offer exponential 
speedups relative to the best classical algorithms known. 
On the other hand, there are results about the limits of 
quantum computation Q , as well as results showing that 
certain problems are hard. For instance, estimating the 
ground state eigenvalue of arbitrary local Hamiltonians 
is a QMA complete problem 

Although there are fundamental problems in complex- 
ity theory that remain open, there is a distinct category 
of problems for which quantum computers can offer sub- 
stantial speedups relative to classical computers. This 
includes problems, such as multivariate integration, path 
integration and multivariate approximation, that suffer 
from the curse of dimensionality in the classical deter- 
ministic worst case. Quantum computers can vanquish 
the curse; see e.g. [5rtZJ- R. E. Bellman introduced 
the term curse of dimensionality referring to multivariate 
problems whose complexity grows exponentially with the 
number of variables and so are impossible to solve when 
the number of variables is large. 

An important problem in physics and chemistry that 
falls in this category is the estimation of the ground state 
eigenvalue of certain multiparticle systems evolving ac- 
cording to the time-independent Schrodinger equation. 
Solving such problems on a classical computer in the 



worst case has cost exponential in the number of parti- 
cles. In particular, the number of state variables d is pro- 
portional to the number of particles and the cost to solve 
the problem with relative accuracy e may grow as e~ d . 
For these reasons researchers have been experimenting 
with quantum computers to solve eigenvalue problems in 
quantum chemistry with very encouraging results 0, 0] ■ 
See also lfl 111. 



In this paper we study a ground state eigenvalue prob- 
lem and we exhibit a quantum algorithm that achieves 
relative error e with cost Cde~^ +S \ where 5 > is an 
arbitrarily small positive number. The cost includes the 
number of queries plus all other quantum operations. 
The algorithm uses C'd log e -1 qubits. The constants C 
and C as well as all constants in our estimates through- 
out this paper are independent of d and e. 

We stress that we are not dealing with an arbitrary 
eigenvalue estimation problem. In our case we are able 
to obtain efficiently a rough but very useful approxima- 
tion of the ground state eigenvector. Abrams and Lloyd 
12| were the first to demonstrate the advantages of ap- 
proximate eigenvectors in solving problems of physical 
interest. Consequently, the cost to implement and simu- 
late the evolution of the Hamiltonian for the amount of 
time prescribed by the accuracy demand determines the 
cost to approximate the ground state eigenvalue. 

We now consider the problem in more detail. If the 
potential is a function of only state variables then the 
ground state energy is given by the smallest eigenvalue 
E\ of the equation 

(-A + V)*i(a;) = #1*1 (a;) for all x e I d := (0, l) d , 
*i(x) = for all x £ dl d , 

with a normalized eigenfunction *i. For simplicity we 
assume that all masses and the normalized Planck con- 
stant are one. The boundary conditions are for particles 
in a box. 

This eigenvalue problem is called the time-independent 
Schrodinger equation in the physics literature and the 
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Sturm-Liouville eigenvalue problem in the mathematics 
literature. We want to approximate E\ with relative er- 
ror e. 

Here, A is the d-dimensional Laplacian and V > is a 
function of d variables. The dimension is proportional to 
the number of particles, e.g. d — 3p. For many applica- 
tions the number of particles p and hence d is huge. We 
consider algorithms that approximate E\ using finitely 
many function evaluations of V. Moreover, we assume 
that V and its first order partial derivatives dV/dxj, 
j = 1, . . . , d, are continuous and uniformly bounded by 1. 

Classical algorithms. Decades of calculating ground 
state eigenvalues of systems with a large number of par- 
ticles have suggested that such problems are hard. We 
sketch a proof that the cost of classical deterministic al- 
gorithms that approximate eigenvalues in the worst case 
grows exponentially with the number of variables. 

Indeed, consider a potential function V and let V be 
a perturbation of V. Then the eigenvalue E\ (V) corre- 
sponding to V and the eigenvalue Ei{V) corresponding 
to V are related according to the formula 



Ex{V) = Ex{V) 



(V(x) - V(a;))#i(a;;V) dx 



0(\\V-V\ 



2 ) 



where ^i(-;V) denotes the eigenfunction corresponding 
to Ei(V). This implies that approximating E\ is at least 
as hard as approximating a multivariate integral in the 
worst case. As a result, any classical deterministic algo- 
rithm for the eigenvalue problem with accuracy e must 
use a number of function evaluations of V that grows as 
e~ d ; see [l3| for details. 

Finite differences are often used for approximating 
Ei. The discretization of the operator — A + V with 
mesh size h = (m + yields an m d x m d matrix 

M h := M h {V) = -Ah + V h . Then one solves the corre- 
sponding matrix eigenvalue problem MhZh.i = Eh,iZh,i- 
Note that A^ denotes the discretization of the Laplacian 
and Vh is a diagonal matrix whose entries are evaluations 
of the potential V at the m d grid points. Mh is symmet- 
ric positive definite and sparse and has been extensively 
studied in the literature 14 , ]j| . For V that has bounded 
first order partial derivatives, using the results of sua 
we conclude 



\Ex-EhA < adh 



(1) 



If Eh t \ is such that \Eh,i — Eh,i\ < C2dh, we have relative 
error 

\1-E h<1 /E 1 \<dh, 

where c' is a constant. The inequality follows by ob- 
serving that Ei is bounded from below by the smallest 
eigenvalue Adh~ 2 sin 2 (7r/i/2) of the discretized Laplacian. 

Quantum algorithm. First we discuss our algorithm in 
general terms and then we provide a complete analysis. 



The key observation is that the procedure we outlined 
above can be implemented on a quantum computer with 
cost that does not grow exponentially with d. This is 
accomplished by modifying quantum phase estimation, 
a well known quantum algorithm for approximating an 
eigenvalue of a unitary matrix W, see e.g., [3, p. 225]. 

Sketch of the algorithm 

1. Consider the discretization Mh — —Ah + Vh of 
— A + V and let h be the largest mesh size lead- 
ing to the desired accuracy, i.e., h < e. The matrix 



W 



JM h /(2d) 



is unitary since Mh is Hermitian. 

2. For W use phase estimation to approximate the 
phase corresponding to e lE h,i/( 2d ) with the follow- 
ing modifications: 

(a) Use the approximate eigenvector 

|0)® 6 |^i)® d 

as an initial state, where is the ground 

state eigenvector of — A^ and can be imple- 
mented efficiently; see the discussion following 
(J4j) below for details. 

(b) Replace W 2 , t = 0, . . . , b — 1, that are re- 
quired in phase estimation, using approxima- 
tions given by high order splitting formulas 
that deal with the exponentials of — A^ and Vh 
separately and can be implemented efficiently; 
see the discussion leading to (|6|) below for de- 
tails. 

The effect of the modifications is to somewhat decrease 
the success probability while increasing the cost of phase 
estimation. Nevertheless, the resulting success probabil- 
ity is at least | , and the cost for implementing the initial 
state and the approximate powers of W does not suffer 
from the curse of dimensionality. 

Theorem 1 Phase estimation with an approximate ini- 
tial state and approximate powers of W with probability 
at least I yields an estimate of Ei with relative error e 



id total cost 



Cde-( 3+5 \ 



for any 5 > 0, using C'dloge 1 qubits, where C and C 
are constants. 

Next we discuss the details of our algorithm and this 
will lead us to the proof of the theorem. The eigenvalue 
of W that corresponds to E h ,x is e lEh - l/{2d) = e 2 ™ Vl , 
where 

<Pi = E htl /(4Trd) 
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is the phase and belongs to the interval [0,1) since E^i < 
Adh~ 2 sin 2 (vrV2) + 1 < dir 2 + 1. 

Quantum phase estimation approximates the phase tpi 
with 6-bit accuracy, where b = \log 2 h~ 1 ~\. The output 
of the algorithm is an index j <E [0, 2 b — 1] such that 
\<Pi - j2~ b | < 2- b . Hence, 



\E h>1 -4:ndj2- b \ < c 2 dh. 
Combining ([T]) and ((2]) we conclude 



Ei - indj 2 b | < c\de + c 2 de = cde. 



(2) 



(3) 



Hence the algorithm approximates the ground state 
eigenvalue E\ by 



E h ,i := 4ndj 2 



-6 



This estimate holds with probability at least (see, e.g., 
[l9j j ) assuming: 

• The initial state of the algorithm is \0)® b \zh,i), 
where \zh,i) is the eigenvector of Mh that corre- 
sponds to Eh i- 



We are given the matrix exponentials W 
0,...,6-l. 



2* 



t = 



In our case, however, we do not know \zh,i) and we 
use an approximation. Similarly, we use approximations 
of the W 2 , t = 0, . . . , b — 1, to simulate the evolution 
of the quantum system that evolves with Hamiltonian 
H = Mh/(2d). We will compute the cost to implement 
these approximations so that ((3]) holds. All these approx- 
imations affect the estimate \ of the success probability 
of phase estimation, but only by a small amount. 

The initial state of our algorithm is 



|o)®Vi 



(4) 



where \ipi)® d is the ground state eigenvector of the dis- 
cretized Laplacian. We know [14[ that the coordinates of 
\ipi) are 



2hsm(jnh), j = l,...,m, 



and \ipi)® d has unit length, 
sented by log 2 m d — 0(dlog 2 



The state is repre- 

ss 1 ) qubits and can be im- 



plemented with d- (3(log 2 s 1 ) quantum operations using 



the Fourier transform; see e.g., [20|, |2l|. We point out 
that here and elsewhere the implied constants in the big- 
O and notation used later in the paper are independent 
of d and e. (From a practical standpoint, it is possible to 
further reduce the cost of the initial state using the algo- 
rithm in [22| but we do not pursue this alternative since 
the analysis of the algorithm becomes more involved.) 

Expanding |V , i)® <i using the eigenvectors of we 
have 



fe=i 



dkUh. 



The approximate initial state reduces the success proba- 
bility of phase estimation by a factor equal to the square 
of the magnitude of the projection of \^p\)® d onto \zh t i), 
to become 1 | 2 ; see, e.g., 121 1221. D ue to the separa- 
tion of the eigenvalues of Mh from [13J we have that 



|di| 2 >l- 



1 



(3tt 2 -2) 2 ' 



This yields that the success probability of phase estima- 
tion with the approximate ground state eigenvector is at 
least 



1 



1 



(3tt 2 - 2f 



> 



(5) 



Now let us turn to the approximation of the matrix 
exponentials. We simulate the evolution of a quantum 
system with Hamiltonian H = Mh/(2d) for time 2*, 
t = 0, 1, 1. Let H = Hi + H 2 where H± = 

—Ah/ (2d) and H 2 = Vh/{2d). Recall that we have chosen 
the largest mesh size such that h < e. The eigenvalues 
and eigenvectors of the discretized Laplacian are known 
and the evolution of a system with Hamiltonian Hi can 
be implemented with d- (9(log 2 e _1 ) quantum operations 
using the Fourier transform in each dimension; see e.g., 
[l8L p. 209]. The evolution of a system with Hamilto- 
nian H 2 can be implemented using two quantum queries 
and phase kickback. The queries are similar to those in 
Grover's algorithm [l8| and return function evaluations 
of V truncated to 0(loge _1 ) bits. 

In particular, we use a splitting formula of order 2fc + l, 
k > 1, to approximate W 2 * = e l{Hl+H ^ 2t by 



Nt 

n< 

£=1 



iAizi 



(6) 



where Ag 6 {H\,H 2 } and suitable zt that depend on t 
and k as described in [23, 24 1. 

From [2^] (see p. 2 using m = 2), the number Nt of 
exponentials needed to approximate H^ 2 by a splitting 
formula of order 2k + 1 with error e i} t = 0, . . . , b — 1, is 



Si Su //, 2' (| 



fc-1 



^2 l \\H 2 \ 



1/(2*) 



for any k > 1. The total number of exponentials required 
for the approximation of all the W 2 is bounded from 
above as follows 



b-l 



N = J2 Nt - 16e ll#il 



t=o 



6-1 

E 2 * 



25 



fc-i 



(8e||ff 2 ||) 



l/(2fc) 



l/(2fc) 



(7) 



25 



< 16e\\Hi\\2° ( - ] (160e2°||ff 2 | 



,1/(2*) 
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where we obtained the last inequality by setting e t 



t = 0, 1. It is easy to check that Ylt=o e * — 

^j. Thus the success probability of phase estimation can 
be reduced by twice this amount [18j, p. 195]. Using (O 
we conclude our algorithm succeeds with probability at 
least 



40 



1 



1 



(3tt 2 - 2) 2 
Since IliTiH < < 2e~ 2 



1 2 
10 - 3' 



— and ||JJ 2 || < l/(2d), the 

algorithm uses a number of exponentials of Hi and H 2 
that satisfies 



N < 32e 



N < C 



-80e\V(2*) /25 
Since we have chosen 2 b = (1/e) we obtain 
d J I 3 



fe-i 



2 6(l+l/(2fe))_ 



for any fc > 0, where C is a constant. 

The optimal k*, i.e., the one minimizing the upper 
bound for JV in (J7), is obtained in [25J, Sec. IV] and 
is given by 



1 , 80e 2 b 

^log 2 5/3— j— 




as de —> 0, 



The number of exponentials corresponding to k* satisfies 



N* 



as de —> 0. 



(8) 



We remark that of the N* matrix exponentials half in- 
volve Hi and the other half involve H2; see the detailed 
definition of the high order splitting formula 23, |24| . 
Since each exponential involving H 2 requires two queries 
the total number of queries is also N* . 

Hence, the number of quantum operations, excluding 
queries, to implement the initial state, the matrix expo- 
nentials involving Hi and the inverse Fourier transform 
yielding the final state of phase estimation is 

0(<f log's -1 ). 



N* 



(9) 

Equations (0), © and flS) yield that the total cost of 
the algorithm, including the number of queries and the 
number of all other quantum operations, is 

Cde- (3+s \ 

where 6 > is arbitrarily small and C is a constant. 
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